#include "op/reduce/reduce.h"

const size_t g_reduceBlockSize = 1024;

__global__ void reduce_sum_double_kernel(const double* input, double* output, int size)
{
    // 首先每个线程内reduce
    int thread_index = threadIdx.x;
    int block_idx = blockIdx.x;
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid >= size) return;

    __shared__ double temp[g_reduceBlockSize];
    temp[thread_index] = input[tid];
    // 同步
    __syncthreads();

    // 线程之间reduce
    size = blockDim.x < size ? blockDim.x : size;
    for (int offset = 1;offset < size; offset *= 2) {
        if (thread_index % (offset * 2) == 0) {
            temp[thread_index] += temp[thread_index + offset];
        }
        __syncthreads();
    }
    
    if (thread_index == 0) {
        // printf("block index: %d, value: %lf\n", block_idx, temp[0]);
        output[block_idx] = temp[0];
    }
    // printf("thread index: %d, value: %lf\n", thread_index, temp[thread_index]);
    
}

void cuda_reduce_sum_double(const double* input, double* output, int size,
    int grid_size, int block_size)
{
    int minGridSize;
    cudaOccupancyMaxPotentialBlockSize(&minGridSize, &block_size, (void*)reduce_sum_double_kernel, 0, size);
    grid_size = (size + block_size - 1) / block_size;
    reduce_sum_double_kernel<<<grid_size, block_size>>>(input, output, size);
}